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A TIME -DEPENDENT METHOD FOR CALCULATING SUPERSONIC 
BLUNT-BODY FLOW FIELDS WITH SHARP CORNERS 
AND EMBEDDED SHOCK WAVES 

By Richard W. Barnwell 
Langley Research Center 

SUMMARY 

A time -dependent numerical method for calculating supersonic blunt -body flow 
fields with sharp corners and embedded shock waves is presented. Axisymmetric and 
symmetric plane bodies can be treated. The method of characteristics is used at the bow 
shock wave and body surface, and a two-step finite -difference method of second-orde* 
accuracy is used between the shock and body. A stability analysis of the finite -difference 
equations which accounts for the effects of a nonorthogonal coordinate system is presented. 

Calculations for flow over bodies with sharp and rounded corners are compared with 
experimental and other theoretical results. The influence of a parabolic free-stream 
velocity distribution on the flow over a spherically blunted cone is investigated, and a 
study of the effect of Mach number on the flow about a truncated cone is made. 

INTRODUCTION 

Time -dependent methods provide a means of calculating supersonic flow past blunt 
bodies with forward-marching numerical techniques because the problem of transient 
flow has hyperbolic equations and is well posed as an initial -boundary -value problem. 
Results for steady flow are obtained from the asymptotic solution to the transient problem. 
These methods can be used to calculate blunt -body flow fields with complicating features 
such as sharp corners on the body contour, embedded shock waves, and large transonic 
regions. 

A number of time -dependent methods have been developed for calculating flow past 
blunt bodies traveling at supersonic speeds. Evans and Harlow (ref. 1) introduced the 
particle -in-cell method, a finite-difference method of first-order accuracy in the mesh 
spacings which treats the bow shock wave as a continuous compression and which is 
formulated in a mixed Eulerian-Lagrangian frame of reference. Eulerian methods of 
first-order accuracy which treat the bow shock wave as a continuous compression have 



been developed by Rich and Blackman (ref. 2); Bohachevsky and Rubin (ref. 3); Gentry, 
Martin, and Daly (ref. 4); Barnwell (ref. 5); Rusanov (ref. 6); and others. 

A finite -difference technique of second-order accuracy which treats shock waves as 
continuous compressions has been developed by Lax and Wendroff (ref. 7). This technique 
has been used by Burstein (ref. 8), Magnus and Gallaher (ref. 9), and others to solve the 
blunt -body problem. 

Another approach to the blunt -body problem is to treat the bow shock wave as a dis- 
continuity. Sauerwein (ref. 10) used this approach and calculated the flow field by using 
the time -dependent method of characteristics. Godunov, Prokopov, and Zabrodin (ref. 11) 
developed a finite -difference method of first-order accuracy which treats the bow shock 
wave as a discontinuity. This method has been extended and used to calculate complicated 
flow fields by Masson, Taylor, and Foster (ref. 12) and McNamara (ref. 13). Other 
methods which treat the bow shock wave as a discontinuity have been developed by Rusanov 
(ref. 14), Xerikos and Anderson (ref. 15), and Moretti and Abbett (ref. 16). The method 
of reference 16 employs the time -dependent method of characteristics at the bow shock 
and body surface, and a finite -difference technique of second-order accuracy between the 
shock and body. 

The purpose of this paper is to extend the method of Moretti and Abbett (ref. 16) so 
that flow fields about bodies with sharp corners and flow fields containing embedded shock 
waves can be calculated. In order to facilitate calculations near sharp corners, a system 
of body -oriented coordinates which focuses at the corners is used, and a conserved form 
of the governing equations with derivatives which are well-behaved near the corners is 
employed. The use of the governing equations in conservation form also permits the cal- 
culation of flow fields containing embedded shock waves. An exact treatment is employed 
to determine the multiple -valued solution at the sharp corner. It should be noted that a 
brief description of the present method and some results calculated with it have been pre- 
sented in reference 17. 


SYMBOLS 

A^Bj^CpD^ quantities defined by equations (5) 
a speed of sound, also radius of flat nose 

B,C matrices defined by equations (A2b) and (A2c) 

D vector defined by equation (A2d) 
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matrix defined by equation (A12) 


eigenvalue of matrix E 
U, U + a, or U - a 
quantities used in equation (B3) 
quantity defined by equation (8) 
matrix defined by equation (A8) 

quantity defined by equation (19), also eigenvalue of matrix G 

total enthalpy, — ^ — (u^ + v^) 

geometry index defined by equation (6) 

curvature of reference line 

curvature of shock 

wavelength of solutions in X- and Y-directions, respectively 

Mach number 

integer 

static pressure 

stagnation pressure 

distance from sharp corner 

radius of curvature of body contour segment 

perpendicular distance from axis of symmetry 


base radius 


radius of curvature of corner 


radius of curvature of nose 
quantities defined by equations (A6) 
arc length of body contour segment 
distance along surface from axis 
time 

velocity components tangent to and normal to shock, respectively 
quantity defined by equation (A14) 

velocity components normal to and tangent to body normals, respectively 

magnitude of velocity 

vector defined by equation (A2a) 

vector used in equation (A5) 

distance along reference line from axis 

normalized distance from body surface defined by equation (9) 

distance normal to reference line 

distance between reference line and body surface 

distance along axis from stagnation point 

exponent used in equation (B3), 0 H a 

complex exponent 


shock angle 



ratio of specific heats 
mesh spacing for time 
mesh spacings for X- and Y -coordinates 
mesh spacings for 77 - and ^-coordinates 
thickness of shock layer 

aa 

9r 
95 

ax 

percent deviation of free-stream velocity at distance r = r^ from axis from 
center -line value 

coordinates along and normal to shock, respectively 
- 7 _ 2tt AY 2 7T AX 

/ 7> S T 9 T 

J-^y -L«X 

6 angle between body normal and free-stream direction 

k damping coefficient which satisfies inequalities (15) 

A quantity defined by equations (A 6 ) 

A x scale factor for X -coordinate 

A^ scale factor for 77 -coordinate 

ix quantity defined by equation (7) 

p density 

a angle between normal to body surface and normal to shock wave 

cp angle defined by equation (A13) 


Y 

At, At 
AX, AY 
Ar/, A | 
6 


• 95 

6 " FT or 


5 ' = 


96 

9x 


or 
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angle defined by equation (A23) 


* 

Subscripts: 

A,D at base points of bicharacteristics 

b at body surface 

C at corner for x = x„ 

c,min 

c at corner 

<t on center line 

l , m indices for grid intersections for X- and Y -coordinates, respectively 

max maximum value 

min minimum value 

s immediately behind shock 

°° in free stream 

Superscripts: 

k time index 

* denotes conditions for smallest value of x where u = a 

iV 

( ) denotes first-order solution 

ANALYSIS 

The present method for calculating numerical solutions for time -dependent, inviscid, 
axisymmetric or symmetric plane flow past blunt bodies traveling at supersonic speeds is 
described in this section. The basic approach is the same as that of Moretti and Abbett 
(ref. 16) in that the method of characteristics is used at the bow shock wave and the body 
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surface, and a second-order finite -difference method is used between the shock and sur- 
face. However, a different coordinate system is used in the present method, and a 
broader class of problems can be treated. 

The problem of transient supersonic flow past blunt bodies is treated as an initial- 
boundary -value problem and is calculated with forward-marching techniques. From an 
initial solution which can be very approximate, the flow is calculated at successive time 
intervals in a closed region which is bounded by the bow shock wave, the body surface, 
the axis of symmetry, and a line extending from the body to the shock which is located 
downstream of the sonic line. At each time step, the flow is calculated first at points on 
the bow shock, then at points on the body, and finally, at points between the body and 
shock. Results for steady flow are obtained after many time steps when the time deriva- 
tives of the flow properties are sufficiently small. 

Calculations at Points Within the Shock Layer 

The basic coordinate system which is used in this paper is curvilinear and is 
oriented with respect to the body as shown in figure 1. The coordinate x is the distance 
along a reference line which is located a constant distance y^ from the body surface, 
and the coordinate y is the perpendicular distance to the reference line and is positive 
outward. The velocity components u and v which are used in this paper are normal 
to and tangent to the body normals, respectively, as shown in figure 1. 

The geometry of the reference line is expressed in terms of its curvature Kj, and 
the angle 6 between the normal to the reference line and a line parallel to the axis of 
symmetry. This angle is shown in figure 1. The quantities K r and 6 are related by 
the equation 


K r = 


de 

dx 


( 1 ) 


The quantity r is the perpendicular distance from a point in the flow field to the axis of 
symmetry, and A x is the scale factor for the x-coordinate. These quantities are given 
by the equations 


r x=x 

r = \ cos <9(x) dx + y sin 9(x) 
J x=0 


( 2 ) 


and 


A x = 1 + yKj. 


( 3 ) 
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The blunt bodies which are treated in this paper have profiles which are constructed 
of segments of constant curvature (straight-line segments and circular arcs). The curva- 
ture of a segment of the reference line is related to the distance y b and the radius of 
curvature of the corresponding body segment by the equation 


1 + 


R 


b 


Straight body segments are associated with parallel straight reference line segments, 
whereas circular arcs on the body are associated with concentric circular arcs on the 
reference line. In the vicinity of a sharp corner, the reference line is a circular arc 
with the radius y^. This arc extends from x c m j n , the value of x on the perpendicular 
to the upstream surface at the corner, to x c max , the value of x on the perpendicular 
to the downstream surface. The length of each segment of the reference line must be 
some integral multiple N of the mesh spacing Ax. The magnitude of Ax is related 
to the number N, the radius of curvature of the body segment R^, and the arc length of 
the body segment S by the equation 



R b + yb 

R b 


If the cross section of the body is circular, the reference line can be constructed with a 
single segment, and the quantities N, S, R^, and y^ can be chosen independently. 

For bodies with two segments with different values of S and R^, the numbers N-^ 
and N 2 of mesh spacings Ax which compose the respective reference line segments 
can be chosen independently, but both sets of quantities S, R^, and N must satisfy the 
equation for Ax. This condition is met if the distance y^ satisfies the equation 


s i S 2 R b,2 + yb 

N 1 R b,l n 2 R b,2 

When the body is constructed of three or more segments, the choice of N, S, and R^ 
for each segment is not arbitrary. 

The governing partial differential equations for inviscid axisymmetric and two- 
dimensional flow can be written in divergence form in terms of the coordinates x and y 
and the time t as 


aAj 9B i 

“at“ + “ex' 


+ 



Di 


0 


(i = 1,2, 3,4) 


(4) 
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In equation (4), the subscript 1 denotes the continuity equation, the subscripts 2 
and 3 denote the tangential and normal momentum equations, respectively, and the 
subscript 4 denotes the energy equation. The quantities A^ B^, Cq, and Dj are 
written as follows: 

A 1 = pX x 

H = P uX x 

A 3 = P vA x 
A 4 = (pH - p)X x 

B 1 = (1 + pj)pu 
Bg = p + (1 + pj)pu 2 
B 3 = (1 + pj)puv 
B 4 = (1 + pj)puH 

> 

C 1 = P vA x 
C 2 = puvA x 

C 3 = (P + P u2 ) A x 

c 4 = pvHA x 

Dj = j(l - p)puA x cos |+ jpvf 

D 2 = j(l - p)pu 2 X x cos |+ puv(K r + jf) 

D 3 = i(l - P)puvX x cos | + jpv 2 f - (p + pu 2 )K r 
D 4 = j(l - p)puHX x cos p+ jpvHf 



where p, p, y, and H are the pressure, density, ratio of specific heats, and total 
enthalpy, respectively. The quantities j, ju, and f depend on the geometry and are 
given by the equations 


and 


J = 


fO (Two-dimensional flow) 
(Axisymmetric flow) 


M = 


(On the axis) 
10 (Off the axis) 


K~ 


f = 


(On the axis) 


X v sin Q 

(Off the axis) 


( 6 ) 

(?) 


( 8 ) 


In order to prevent the bow shock from cutting across the grid lines, the 
y-coordinate is normalized with the local shock layer thickness 6 as shown in figure 2. 
The normalized coordinate Y, which is given by the equation 


Y = 


y + yb 

6 


( 9 ) 


has the values of zero at the body surface and 1 at the shock wave. The partial deriva- 
tives with respect to the old independent variables t, x, and y are related to the 
partial derivatives with respect to the new variables t = t, X = x, and Y = y by the 
expressions: 

_9__ d_ Y5 9 
9t 9t " 6 9Y 

d_ _ _9_ Y6|_ _9_ 

9x 9X ' 6 9Y 


9 _ 1 9 
9y 5 9Y 

where 6 and 6* are the derivatives of 6 with respect to r and X, respectively. 
Equation (4) can be written in terms of t, X, and Y as 
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(i = 1,2, 3, 4) (10) 



There are several advantages to using the governing equations in this particular 
conservation form. In the first place, it improves the results at grid points adjacent to 
sharp corners because the quantities Cj - Y are proportional to Y so 

that the derivatives of these quantities with respect to Y in equation (10) are well- 
behaved at these points although the individual properties have large Y-derivatives there. 
This condition is true because the reference line subtending the corner is a circular arc 

with the curvature K„ = — so that the scale factor satisfies the equation 

r y b 



and the following relationship holds : 
q -y^AjBH- Bjfi') ccY 

In the second place, the conservation form permits the calculation of flow fields con- 
taining embedded shock waves, no special consideration being given to computations in 
the vicinity of those waves. 

A two-step Lax-Wendroff finite -difference method is used to obtain solutions of 
second-order accuracy at points between the shock wave and body surface. The first 
step of the calculation is the determination of a solution of first-order accuracy in the 
mesh spacing At, Ax, and AY at the new time r from the known solution at time 
t - At, whereas the second step is the correction of the solution at time t to second- 
order accuracy. Let k, l, and m be the indices for the time r and the coordinates 
X and Y, respectively. The first-order solution at time t is obtained with the 
equations 





k-l 


m+1 


k-l 

C^i) + 

V n 7,m-l 


(M 


k-l 


Z+l,m 


+ 


(M 


k-l 

Z-l,m 


+ At 



(i = 1,2,3, 4) 


(ID 


where values for the partial derivatives 9A^/8 t at the point (k - 1, l, m) are determined 
by evaluating the right-hand side of equation (10) at this point. The partial derivatives 
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with respect to X and Y in equations (10) are approximated with central difference 
formulas. Once the quantities Aj are known, the flow properties p, p, u, and v 
can be determined. It should be noted that the scale factor A x , and hence the shock 
layer thickness 6 at time r, is needed in order to determine the density from the 
known value of A^. It is for this reason that for each new time step, the solutions at 
points on the bow shock wave are determined before the solutions at points between the 
shock and body. 

The second-order solution at time r is obtained with the equations 



where the values for the partial derivatives at the point (k, l, m) are deter- 

mined by evaluating the right-hand side of equation (10) at this point with the first-order 
solution just obtained. The terms of fourth order in equations (12) are not physical and 
are added to eliminate instabilities. 

The von Neumann conditions provide a means of estimating the upper bound of the 
time step At which can be used with a given set of finite-difference equations. It is 
shown in appendix A that, in terms of the velocity components u and v and the non- 
orthogonal coordinates X and Y used in this paper, this condition can be written as 


At S 


\/(*x *X) 

- I \ V ^ \ v a W J 

2 + (6 AY) 2 / 

(Y6’AX) 2 1 

t , (Y6* AX) 2 

2 

2(A X AX) (6 AY) " 

| * <» iY > 2 | 

(\ x AX) 2 + (6 AY) 2 

(A x AX) 2 + (6 AY) 2 


(13) 


where a is the local speed of sound. In practice, it has been found that this inequality 
can be replaced by a simpler expression of the form 


At = 


fy + 1 V 


max, 00 


min(A x AX, 6 AY) 


(14) 
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where V max>00 = «/2H^ is the maximum speed in the free stream. It is also shown in 
appendix A that the damping coefficient k in equation (12) must satisfy the inequalities 

0 < * = W ( 15 ) 


Calculations at Points on Bow Shock Wave 

As in reference 16, a time -dependent method of characteristics is employed at the 
shock wave and body surface. The present method differs from that of reference 16 in 
that the characteristic compatibility relation is integrated along a bicharacteristic rather 
than along a quasi -one -dimensional characteristic. A general discussion of the time- 
dependent method of characteristics for two-dimensional flow can be found in 
reference 10. 

An orthogonal curvilinear coordinate system 77 , £ is established at the shock wave 
as shown in figure 3, where 77 and £ are the distances along and normal to the shock, 
respectively, at time r. The velocity components in the 77 - and ^-directions are U 
and V, respectively. The angle a between the normals to the shock wave and the 
reference line satisfies the relation 

cl 

tan a = — (16) 

A x 


The basic problem is to determine the solution for the flow properties and the 
shock speed at points on the shock at time r from the known solution at r - At. The 
equations which are solved at points on the shock wave are the Rankine-Hugoniot rela- 
tions and one characteristic compatibility relation. If the subscripts 00 and s denote 
the quantities in the free stream and immediately behind the shock, respectively, the 
Rankine-Hugoniot relations can be written as 


P s (v s - sec a 6) = - sec a 6) 

P s + P s ( V s - sec a 6) 2 = Poo + Poo (V M - sec a 6) 2 

— r — + ~(V - sec a 6)^ = — — + i(v - sec cr 6) 

y - 1 P 0 2' s ) y-lp„2V°° / 




U S = U oo 


(17) 
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The characteristic compatibility relation is written as 


pa^ = + Pj)pa 2 |^+ j(l - p)pa 2 U ^ cos p 

+ pa 2 V(K s + jg) - KgpaU 2 ] (18) 

where d/dr is the substantial derivative, j and p are given by equations (6) 
and (7), respectively, and g is given by the equation 


g = 



sin (3/r 


(On the axis) 
(Off the axis) 


(19) 


The quantities /3, K s , and A^ in equation (18) are the angle between the normal to the 
shock wave and the free-stream direction, the curvature of the shock, and the scale 
factor for the t? - coordinate, respectively, and are given by the equations 


p = 0 - a 


K 


s 


977 


A^J - 1 + Kg£ 


The compatibility relation (18) is integrated from a point between the shock and body at 
time t - At to the point on the shock at time r along the bicharacteristic which satis- 
fies the differential equations 


d?? _ U 
dr Afj 

££= v + a 
dr 




j 


( 20 ) 


The equation thus obtained gives a linear relationship between the pressure p g and the 
velocity component normal to the shock V . This equation and the Rankine -Hugoniot 
relations (17) are solved simultaneously to determine values for the flow properties p g , 
p g , U g , and V g and the component of the shock velocity normal to the reference line 5 
at time r. 
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The following procedure is used to obtain solutions at the bow shock. First, a 
tentative value for 6(77, t), the location of the shock wave at time t, is obtained from 
the equation 


6(tj,t) = 6(17, t - At) + At 6(77, t - At) 


at all points on the shock wave, the derivative 6' (77, t) is determined geometrically, and 
the angle a is evaluated with equation (16). Then it is assumed that 6(77, t) = 6(77, t - At), 
and an initial solution for p g , p g , U s , and V g at time t is obtained from equa- 
tions (17). 

The point A where the bicharacteristic through the shock point s intersects 
the plane t - At is determined by iteration. In order to accomplish this iteration, 
equations (20) are approximated by the finite -difference expressions 


A 77 _ _1 
At 2 




If = + a >s + < v + a > A ] 






( 21 ) 


Initially, it is assumed that the properties at A have the same values as those at s. 
These values are substituted into equations (21) and a first estimate for the location 
of A is obtained. Improved values for U^, V^, and a^ are determined by inter - 
polating the known solution at time t - At and are used in equations (21) to obtain a 
second estimate of the location of point A. This process is repeated until the results 
for the location of point A converge. 

When point A has been located, the compatibility relation (18) is integrated from 
A to s along the bicharacteristic. In order to perform this integration, the right- 
hand side of equation (18) and the coefficient pa on the left-hand side are evaluated at 
s and A, and the results are averaged. The integrated compatibility equation and the 
Rankine -Hugoniot relations (17) are solved simultaneously for improved values of the 
quantities p g , p g , U g , V g , and 6 at time t, and a second approximation for the 
shock location is determined with the equation 

6(77, t) = 6(77, t - At) + | AtJ6(t7,t - At) + 6(77, t)J 

Once the improved solution has been obtained at all points on the shock wave, new esti- 
mates of the derivative 6'(t7,t) and the shock angle a are determined for each point. 
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In general, the second set of values for the flow properties and shock velocity and 
location at time r differ slightly from the first set. Converged results are obtained 
by repeating the process described several times, beginning with the location of the 
point A, and continuing until the next set of improved values for the flow properties and 
shock wave parameters have been determined. 

Calculations at Points on the Body Surface 

The basic x,y coordinate system shown in figure 1 is used in connection with the 
time t at points on the body surface. Therefore, the governing partial differential 
equations are given by equations (4). These equations are combined linearly to obtain 
characteristic compatibility relations which are used in conjunction with the boundary 
condition 


v = 0 

to obtain solutions for the flow properties at the body surface. The pressure can be 
determined with the relation 


. oa = - i_ 

P dt A x 


dt 


(1 + 4 j)pa 2 ^+ j(l - p)pa 2 u cos 9 + pa 2 v(K r + jf) - K r pau 2 


( 22 ) 


where j, p, and f are given by equations (6), (7), and (8), respectively. Equa- 
tion (22) is integrated along the bicharacteristic which satisfies the differential 
relations 


dx _ u 
dt A x 



j 


(23) 


The density is determined with the compatibility solution 


dP _ a 2 
dt dt 


0 


(24) 


which is integrated along the streamline which wets the body surface. This streamline 
has the slope 


dx _ u 
dt A x 


( 25 ) 
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in the x,t surface. The surface velocity is determined either from the transient 
Bernoulli equation 


dH 1 9£ = o 
dt ' p 9t 


( 26 ) 


which is integrated along the surface streamline or by the compatibility relation 


pa*. -,»*(£♦ J.SO*) 


(27) 


which is integrated along the bicharacteristic in the x,t surface with the slope 


§ = Sf* ( 28 ) 

at 

An iterative procedure is used to determine the solution at the body surface at 
time t which is similar to that employed at the shock wave. For each cycle of the 
iteration, values for the flow properties are determined at all the points on the body sur- 
face. The iteration is continued until the solutions for successive cycles of the iteration 
converge. In order to start the first cycle, the values of the flow properties at time t 
are set equal to the known values at time t - At. 

The basic procedure which is employed at the body surface consists of using equa- 
tions (22), (24), and (26) or (27) to determine the pressure p, the density p, and the 
tangential velocity component u, respectively. This procedure is used to determine the 
solution for all values of x for which the surface is continuous and for x = Xg^j^ 
if the flow approaching the corner is supersonic. 

The location of point A, where the bicharacteristic through the body point b at 
time t satisfies the differential relationships (23) and intersects the surface t - At, 
is determined iteratively in the same manner that was used at the shock wave. The 
compatibility relation (22) is integrated from A to b to determine the pressure p^. 
Then the intersection point B of the surface t - At and the streamline bicharacteristic 
through b, which satisfies equation (25), is determined and equation (24) is integrated 
from B to b to obtain the density p^. At this stage either equation (26) or equa- 
tion (27) is integrated to determine the tangential velocity u^. It should be noted that at 
the axis of symmetry, the condition % = 0 is imposed. When the transient Bernoulli 
equation (26) is used, the partial derivative 8p/9t is evaluated explicitly. If the x- 
and t-coordinates of the points D and b are denoted by x^, t - At, and x^, t, 
respectively, the finite -difference expression for this derivative is written as 
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8p = 1 P(x b> t ) - P(x b> t - At) p(x D ,t) - p(x p > t - At) 

9t 2 At At 

The density p in equation (26) is approximated with the average of the densities at the 
points D and b. 

If the flow approaching a sharp corner on the body contour is subsonic, an alternate 
procedure is used to determine the solution at the corner for x = x C)inin . This proce- 
dure consists of simultaneously solving the integrated energy equation (24) and either the 
integrated transient Bernoulli equation (26) or the integrated compatibility equation (27) 
subject to the sonic condition u = a. 

It is shown in appendix B that the solution at the corner for 

v . < y < v 

A c,mm - A - A c,max 

is the transient analog to the Prandtl -Meyer solution, which is given by equations (B5) 
and (B6). 


Initial and Boundary Conditions 

A complete initial solution can be constructed from a specified shock shape and 
surface pressure distribution, both of which can be very approximate. In this paper, the 
initial shock shape is a conic section, and the surface pressure distribution for smooth 
bodies is the Newtonian distribution. If the body has a sharp corner where the flow is 
sonic, an empirical distribution with a singularity at the corner is used. It is assumed 
that the initial shock velocity is zero so that the flow properties at the shock wave can be 
determined completely with the Rankine-Hugoniot relations and the approximate shock 
shape. The density and tangential velocity at the surface are determined with the approx- 
imate surface pressure and the principles of conservation of surface entropy and total 
enthalpy. The initial properties at points between the body and shock are obtained by 
interpolating linearly between the body and shock along normals to the body surface. 

As discussed previously, the boundary conditions applied at the shock and body are 
the Rankine-Hugoniot relations for a moving shock and the flow -tangency condition, 
respectively. The boundary condition at the axis of symmetry is that the tangential 
velocity vanish there. The flow properties at the downstream boundary are determined 
by extrapolation. If this boundary is located in a region where the flow is supersonic, the 
errors incurred by the extrapolation do not propagate back upstream and hence do not 
influence the calculation appreciably. 
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RESULTS AND DISCUSSION 


All the results which are presented are for a perfect gas with a ratio of specific 

heats of 7/5 and for axisymmetric bodies at zero angle of attack. Unless otherwise 

specified, the free stream is uniform. When it is feasible, the pressure distributions 

_o 

are normalized with the quantity instead of the stagnation pressure so that the 

calculated value of the stagnation pressure can be presented. Solutions are presented 
in detail for several cases and are compared with experimental and other theoretical 
results in order to demonstrate the accuracy of the method. In addition, results are 
presented for flow fields with nonuniform distributions of the free-stream properties 
and flow fields with overexpansion and recompression regions. 

Sample Calculations 

Results are presented for the flow fields about two flat -face cylinders, one with a 
rounded corner and the other with a sharp corner, and a blunted cone with a sharp corner. 
For each of the flat -face cylinders, detailed comparisons are made between experimental 
results and the result of two calculations with different mesh spacings. In the case of the 
blunted cone, comparisons are made between experimental results and the results of the 
present and other methods. 

Flat -face cylinder with rounded corner .- The ratio of the radius of the corner to 
the base radius is 0.4, and the free-stream Mach number is 2.49. The shock-wave and 
sonic -line locations are shown in figure 4, the surface pressures and velocity distribu- 
tions are given in figure 5, and flow property profiles along several normals to the body 
surface are presented in figure 6. The results for solution 1 of the present method 
depicted in figures 4, 5, and 6 were calculated with a grid with 21 mesh spacings in the 
X-direction and 6 in the Y-direction. The results for solution 2 were calculated with a 
grid with 14 and 4 mesh spacings, respectively, in the X- and Y-directions. 

The experimental data for the shock -wave location shown in figure 4 were obtained 
coincident with the data presented in reference 18. The experimental data for the pres- 
sure and velocity were obtained from reference 18 and are depicted in figure 5 with a 
solid line as a function of the distance s along the surface from the axis. The normal 
to the surface designated as I in figure 6 is coincident with the stagnation streamline, the 
normal designated as II intersects the body on the shoulder and is inclined at an angle 
of 30 ° to the free-stream direction, and the one designated as in intersects the body at 
the junction of the shoulder and the afterbody. 

Flat -face cylinder with sharp corner .- The Mach number for this case is 2.81. The 
shock -wave and sonic -line shapes are shown in figure 7, the surface pressure and velocity 
distributions are given in figure 8, and the flow property profiles along the stagnation 


19 



streamline and the lines normal to the upstream and downstream portions of the body 
surface at the sharp corner are presented in figure 9. The experimental data shown in 
figures 7 and 8 were obtained from results given in reference 19. In figure 9, the normal 
to the body along the stagnation streamline is designated as I; the normal to the upstream 
face at the corner, where X = X c min , is designated as II; and the normal to the down- 
stream face at the corner, where X = X c>max , is designated as III. 

It is seen in figures 7, 8, and 9 that the numerical solutions differ slightly as a 
result of different mesh spacing sizes. For example, the numerical values for the pres- 
sure on the face near the corner are higher than the experimental values, whereas the 
values for the velocity on the face are generally lower than the experimental results. 

Both of these trends decrease as the size of the mesh spacings is decreased. However, 
it can be seen in figures 8 and 9 that the two solutions have the same singular character 
in the vicinity of the sharp corner. Note that at a sharp corner, the velocity component v 
shown in figure 9(d) is zero only where x = x c>m ^ n . (See appendix B.) 

Blunted cone with sharp corner.- The half-angle of the cone shown in figure 10 is 
60°, and the nose of the cone is spherically blunted and the ratio of nose radius to base 
radius is 0.25. 

In figure 10(a), the shape of the bow shock wave, as predicted by the present method 
for a free-stream Mach number of 4.63, is compared with that predicted experimental data 
(ref. 20) and by the one-strip method of integral relations (ref. 21). It is seen that the 
results of the method of integral relations are in better agreement with experiment than 
the present results in the stagnation region. In the vicinity of the sharp corner, both 
methods tend to underestimate the thickness of the shock layer. The results for the sur- 
face pressure distribution are compared with experiment in figure 10(b). As noted in ref- 
erence 21, the method of integral relations tends to underestimate the pressure on the 
conical part of the body. It is seen that the present results show good agreement with 
experiment in this region. 

The results of the present method and the method of integral relations for the shock 
shape and surface pressure distribution for a free-stream Mach number of 10 are com- 
pared in figure 11. Also shown in the figure are the results of Masson, Taylor, and 
Foster (ref. 12) for a free-stream Mach number of 9. As has been mentioned previously, 
the results of reference 12 were calculated with the Godunov method, which is a time- 
dependent method of first-order accuracy. Because results for different Mach numbers 
are being compared in figure 11, the pressure distributions are normalized with the 
stagnation -point pressure. It is seen that the present method predicts that the shock- 
layer thickness is less than that predicted by the other two methods at all points along 
the body and that the surface pressure is greater than that predicted by the other methods. 
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The differences in the free -stream Mach numbers for the results shown in this figure do 
not account for the differences in the results since the Mach numbers are so large. In 
fact, the pressure predicted by the Godunov method on the conical part of the body for a 
free-stream Mach number of 9 is less than that predicted by the present method for a 
free -stream Mach number of either 4.63 or 10. 

It is seen in figures 10(b) and 11(b) that the present results for the surface pressure 
are irregular just upstream of the corner. This irregularity decreases as the grid is 
refined and is probably related to the mesh-spacing size effect observed for the flat -face 
cylinder with the sharp corner. 

Calculations for Nonuniform Free Stream 

The magnitude of the free-stream velocity is assumed to vary parabolically with 
the perpendicular distance from the axis so that flow fields about bodies in a nonuniform 
free stream such as that found in some wind tunnels or in the wake of a forebody can be 
investigated. The expression for the velocity is written in the form 

V„(r)=V„ t |i- £ (r/r b ) 2 J 

where ^ is the magnitude of the free-stream velocity at the center line, r^ is the 
base radius of the body, and e is the percent deviation of the free-stream velocity at a 
distance r = r^ from the axis from the value It is assumed that the free-stream 

pressure and total enthalpy are constant and that the flow in the free stream is parallel 
to the axis. This type of nonuniformity leads to a strong dependence of the dynamic pres- 
sure on the distance from the center line and becomes more pronounced as the center-line 
Mach number is increased. 

The body which is treated is a spherically blunted cone with a sharp corner, a half- 
angle of 60°, and a nose radius to base radius ratio of 0.25. The Mach number at the 
center line is 10. 

Bow shock wave and sonic line shapes .- The locations of the shock waves and the 
sonic lines resulting from three nonuniform free streams are compared with those for 
the uniform free stream in figure 12. It can be seen that the nonuniformity has a marked 
influence everywhere in the flow field. The general trend as e is increased is for the 
shock wave to move closer to the body in the stagnation region and further from the body 
in the transonic and supersonic regions. The sonic point at the bow shock moves outward 
away from the axis as e is increased. The results indicate that the sonic point on the 
surface is positioned at the corner for e = 0, 0.005, and 0.01, but that it is located 
upstream of the corner for e = 0.03. This displacement of the sonic point will be dis- 
cussed subsequently. 
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Also shown in the figure are the results of the one -strip method of integral rela- 
tions for e = 0.01. The results were calculated with a modified form of the method 
described in reference 21. 

Surface pressure dis tribution. - In figure 13, the surface pressure, which is made 
—2 

nondimensional with is plotted as a function of the distance along the surface 

from the axis. It can be seen that for values of e of 0.005, 0.01, and 0.03, the pressure 
near the corner (s/r^ ~ 0.8) is reduced by about 10, 20, and 50 percent, respectively. It 
should be noted that these pressure reductions are accompanied by equally severe reduc- 
tions in the free -stream Mach number. At a distance r = r^ from the axis, for 
e = 0.005, the free-stream Mach number is reduced by 9 percent from 10 to 9.1; for 
6 = 0.01, it is reduced by 16 percent to 8.4; and for e = 0.03, it is reduced by 34 percent 
to 6.6. 

The results of the one -strip method of integral relations for e = 0.01 are also 
shown in the figure. The close agreement of the results of the present method and those 
of the method of integral relations is encouraging in the absence of experimental data on 
the subject. 

The Newtonian pressure distributions are presented in figure 13 for e = 0, 0.01, 
and 0.03. The standard Newtonian equation was used, but the radial dependence of the 

-2 

quantity was included. The results of the present method for e = 0 and 0.01 

differ considerably from the Newtonian results as expected since the flow is subsonic all 
along the face and is singular at the corner. The Newtonian results for e = 0.03 indicate 
that the sonic point is located on the conical portion of the body upstream of the corner at 
approximately the same location as the present results. The two sets of results are in 
close agreement in the supersonic region. 

In general, there is a reluctance to accept solutions which indicate that the sonic 
point on the surface of a body in a symmetric inviscid flow field is located on a straight 
segment of the surface. In fact, Nikolskii and Taganov (ref. 22) have proved that for 
plane potential flow such a solution cannot exist and Shifrin (ref. 23) has extended the 
proof to cover symmetric, nonisentropic plane flow. However, it should be noted that 
this proof does not apply directly to axisymmetric flow. 

It is possible that the truncation error of the present method is responsible for the 
displacement of the sonic point upstream of the corner for the case e = 0.03 since this 
error is viscous-like. However, it is also possible that the present results represent 
the inviscid flow correctly, and that the displacement occurs because the dynamic pres- 
sure at points off the axis for e = 0.03 becomes too small to support subsonic flow at 
the surface all the way to the corner. 
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Calculations for Overexpansion-Recompression Flow Fields 

The bodies about which flow is calculated in this section are truncated cones with 
rounded shoulders between the face and conical afterbody. The flow overexpands around 
this shoulder and recompresses because of the presence of the afterbody. The recom- 
pression starts as a continuous compression on the surface of the afterbody just down- 
stream of the shoulder. Experiments (refs. 19 and 24) show that if the flow has expanded 
to the supersonic state, the Mach lines in the compression region may converge so that 
an embedded shock wave is formed between the surface of the afterbody and the bow shock 
wave. 


Comparisons with experiment and another theoretical method.- The free-stream 
Mach number is 1.79, and the configuration about which the flow is calculated is body 2 
of reference 23. The half -angle of the conical afterbody of this configuration is 30°, and 
the ratios a/r b and r c /r b are 0.344 and 0.2, respectively. The results of the present 
method, the one-strip method of integral relations by Traugott (ref. 25), and the experi- 
mental data of Hastings, Persh, and Redman (ref. 24) for the shockwave and sonic line 
locations, and the surface pressure distribution are shown in figures 14 and 15. The 
pressure distribution in figure 15 is normalized with the stagnation value because an 
experimental value for the stagnation pressure is not given in reference 24. 

The shadowgraph of this flow field indicates that there is a weak embedded shock 
wave in the recompression region over the afterbody. The present results for the dis- 
tributions of pressure in this region along the lines Y = Constant are shown in figure 16. 
The abscissa is the distance from the upstream edge of the conical afterbody. These 
results indicate that the flow is being compressed, but the relative magnitude of the pres- 
sure increase is small and decreases as the distance from the surface is increased. The 
points where the experimentally observed shock wave crosses the lines Y = Constant, 
which are plotted on the figure, are located where the pressure gradient is relatively 
large. 

Effect of Mach number on shock-wav e and sonic -line locations .- The influence of the 
free-stream Mach number on the structure of the shock layer about a 45° truncated cone 
with rounded shoulders is shown in figure 17. The relative magnitudes of the radius of 
the flat face, the corner radius, the length of the conical surface, and the base radius are 
1, 0.25, 2, and 1.664, respectively. The cone is followed by a cylindrical afterbody. 
Results are presented for Mach numbers of °°, 2.8, 2.6, 2.4, and 2.0. 

Also shown in the figure are results for the flat -face cylinder with the same value 
of the ratio r c /a as the truncated cone -cylinder. Comparisons of the results for the 
two bodies show the influence of the conical portion of the truncated cone-cylinder on the 
flow field. 
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For a free -stream Mach number of infinity, the sonic line and bow shock wave in 
the nose regions of the two bodies coincide. Thus, the flow field ahead of the first 
shoulder of the truncated cone-cylinder is not influenced by the presence of the conical 
portion of the body. The flow expands around the first shoulder and becomes supersonic. 
It then compresses on the conical part of the body and becomes subsonic. At the second 
shoulder, it again expands to the supersonic state. Qualitatively, the structure of the 
shock layer for a free -stream Mach number of 2.8 is essentially the same as that for 
M m = oo. However, when M^ = 2.8, the relative size of the subsonic region over the 
conical part of the body is much larger. 

The regions of subsonic flow are connected for a Mach number of 2.6; thus, the 
flow in the stagnation region is influenced by the presence of the afterbody. However, 
this influence is so slight that the locations of the bow shock wave in the nose regions of 
the truncated cone and flat face cylinder are essentially the same. It should be noted 
that a bubble of supersonic flow remains at the first shoulder. 

When the Mach number is 2.4, the zone of subsonic flow over the conical portion of 
the body enlarges so that it extends all the way to the bow shock wave. There are two 
isolated regions of supersonic flow at this Mach number: one at the bow shock wave above 
the first shoulder, and one on the first shoulder. The conical part of the body still has 
little influence on the flow in the stagnation region; thus, the shock-detachment distances 
for the truncated cone and the flat cylinder are the same. 

For a Mach number of 2, the flow field upstream of the second shoulder is subsonic 
with the exception of the supersonic bubble at the first shoulder. It is seen that at this 
Mach number the entire subsonic region is affected noticeably by the presence of the 
conical surface. 

It should be noted that for y = 1.4, the weak shock solution for flow over a pointed 
45° cone with the shock attached at the apex predicts that the flow at all points in the 
shock layer is supersonic when the Mach number is infinity and 2.8, and subsonic when 
the Mach number is 2.6 and 2.4. No attached-shock solution exists for a 45° cone when 
y = 1.4 and M m = 2. 

The compression which occurs over the conical part of the truncated cone starts on 
a continuous compression at the first shoulder. At higher free-stream Mach numbers, 
the compression fan probably merges to become a weak embedded shock wave which 
intersects the bow shock wave. On the basis of the present results, it cannot be deter- 
mined whether the embedded shock forms since the present method smears weak shocks 
so that they cannot be distinguished from continuous compressions. It is clear that the 
continuous compression extends all the way to the bow shock wave for = 2 because 
the flow is subsonic throughout the recompression region. 
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Effect of Mach number on surface pre ssure distribution . - In figure 18, the surface 

pressure distributions, which have been nondimensionalized with the stagnation pressure, 

are plotted as a function of distance along the surface from the axis for the five values 

— 2 

of M^. The stagnation -point pressure p^ rather than the quantity is used for 

normalization because it leads to a better overall correlation of the results. The sonic 
pressure level and the values of the surface pressure for a pointed 45° cone with the 
shock attached at the apex are shown also. As noted before, no attached- shock solution 
exists when = 2. 

It is seen that the results correlate very well in the stagnation region. On the first 
shoulder a Mach number effect appears which becomes very pronounced on the conical 
part of the body. This effect is monotonic, the higher nondimensional pressures corre- 
sponding to the lower Mach numbers. Note that the distribution for M^ = 2 appears to 
provide an upper bound for the other distributions on the rear part of the conical surface, 
and that the other distributions tend to approach the surface pressure value for the 
pointed 45° cone at the appropriate Mach number. Mach number effects are observed on 
the downstream part of the second shoulder but they are small. 

CONCLUDING REMARKS 

A time -dependent numerical method is presented for calculating supersonic blunt - 
body flow fields with sharp corners and embedded shock waves. Both axisymmetric and 
symmetric plane bodies can be treated. The method of characteristics for transient flow 
was used at the bow shock wave and body surface, and a two-step finite -difference method 
of second-order accuracy which employs a shock and body-oriented coordinate system is 
used between the bow shock and body surface. Exact closed form expressions are used 
to determine the multiple -valued solutions at the sharp corners. A stability analysis of 
the finite -difference equations which accounts for the effect of the nonorthogonal coordi- 
nate system is presented. 

Some results for bodies with both sharp and rounded corners were compared with 
experimental data and other theoretical results to establish the accuracy of the method. 

The method was used to investigate the effect of a parabolic free-stream velocity 
distribution on the flow field about a 60° blunted cone with a sharp corner and a ratio of 
nose radius to base radius of 0.25. It was found that this nonuniformity tended to decrease 
the thickness of the shock layer adjacent to the blunted nose and increase the thickness 
adjacent to the conical surface. In addition, it was shown that increasing the nonuniformity 
caused a large decrease in the pressure on the conical part of the body. 

A study of the effect of free-stream Mach number on the flow about a 45° truncated 
cone-cylinder with rounded corners was also made. It was found that for free-stream 
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Mach numbers between 2 and infinity, the flow overexpanded to the supersonic state on 
the corner between the face and the conical surface and then recompressed on the conical 
surface. For Mach numbers of 2.4 and higher, the pressure on the conical surface 
approached the value for similar flow over a pointed cone. It should be noted that similar 
flow does not exist for a pointed 45° cone for free-stream Mach numbers much below 2.4. 
It was found that the pressure distribution on the conical afterbody for M m = 2, when 
normalized with the stagnation value, serves as an upper bound for the normalized distri- 
bution for the larger free-stream Mach numbers. For free-stream Mach numbers of 2.4 
and above, the bow -shock -wave location in the stagnation region is not influenced appre- 
ciably by the presence of the coniGal part of the body. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 9, 1970. 
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APPENDIX A 


STABILITY ANALYSIS FOR FINITE -DIFFERENCE EQUATIONS 

In this appendix, the von Neumann condition for the linearized form of the present 
difference equations is determined. This condition specifies an upper bound for the 
mesh spacing At above which it may be expected that the magnitude of a given infini- 
tesimal error will increase at successive time steps. In general, difference schemes 
which amplify small errors produce solutions which either diverge or contain large 
errors. 

The basic procedure used in this paper to determine the von Neumann condition is 
given by Richtmyer (ref. 26). A technique employed by Van Leer (ref. 27) is used to 
account for the difference in the mesh spacing sizes for the X- and Y -coordinates. A 
nonphysical dissipation function of fourth order similar to that used by Richtmyer and 
Morton (ref. 28) is employed to avoid neutral stability at points where the components of 
velocity vanish or are sonic. This treatment accounts for the nonorthogonal nature of 
the X,Y coordinate system. 

The differential equations (10) can be written in matrix form as 


9w B 8w C 9w p 
8r A x dX 6 8Y 


0 


where the vectors w and D and the matrices B and C 


are written as 
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(A2c) 


(A2d) 


The expanded forms of the finite -difference equations (11) and (12) are 


~k 1 / k-1 

V = 4l w Wr +w ’ 1 - +Wi ~' 1+W 


k-1 .J’ 1 AW k " ] 

m • r W l-l,m +w l,m+l +vf l, 
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(A3) 
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and 


w 


= w* 


-1 1 At 


l , m l,m 4 x x AX 




1 At r (~ k ~k ^...k-l „,k-l \ A 

4 6 AY C \ W I,m+l ” w i,m-l + w f,m+l " w I,m-lJ " A 

- K [( w fc2,m - 4w |Vl,m + H'm ' + w f-2,m) 

+ Kmt2 - 4w t"m + l + 6w i,m ' 4w f,m-l + "f/nU)] 

respectively. In this appendix it shall be assumed that the matrices B and C are 
evaluated at the point (k - 1 , l, m). 

Let the small errors be represented by a series of terms of the form 


(A4) 


w 


= w Q g k exp [i(fi + m 77 )] 


(A5) 


where 


g = exp(o! At) 

l = p- AX 
L x 


V = 


p- AY 


The quantities L x and L y are the wavelengths of the error solutions in the X- and 
Y-directions, respectively, and define the quantities A, r^, and rg as 
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The solution at time r in terms of that at r - At for a small time step At is 
obtained by substituting equation (A5) into equations (A3) and (A4) and taking the limit 
as At approaches zero so that the quantities At D in equations (A3) and (A4) vanish. 
This solution is written as 

w(X,Y,t) = Gw(X,Y,t - At) (A7) 

where the matrix G is given by the equation 
G = ^1 - 4kJ(1 - cos i) 2 + (1 - cos 

- i A 2 ^r^B sin | + rgC sin fj^ 

- i i(2 + cos 1 + cos fj) A^r^B sin | + rgC sin rj^ (A8) 

Equation (A8) is the amplification matrix for the difference scheme given by equations (A3) 
and (A4). 

The solutions at times t and t - At are related also by the equation 


w(X,Y,t) = e a At w(X,Y,t - At) = gw(X,Y ,t - At) (A9) 

It can be seen from equation (A9) that the magnitude of the solution at time t is less 
than that at time t - At; hence, the solution satisfies the von Neumann condition for 
stability if the magnitude of g is less than one. 

It can be shown by combining equations (A7) and (A9) that g represents the eigen- 
values of the matrix G: 


(G - Ig)w = 0 (A10) 

The quantity I is the identity matrix. An expression for g can be obtained from the 
secular equation for equation (A10): 

|G - Ig| = 0 (All) 

Let the matrix E and the quantities cos cp and U be defined as 

E = r^B sin \ + v^C sin r\ (A12) 
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cos cp = 


rj sin 


| - r 2 sin 77^Y6'/A x ] 


y jr^ sin | - r 2 sin ? 7 ^YS'/A x jj 2 + r 2 sin 2 ?) 


(M3) 


and 


U = u cos cp + (v - Y6)sin cp 
Therefore, the matrix E can be written as 


(A14) 
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= /[ r 1 sin « " r 2 sin ^(Y^Ax)] 2 + r 2 



U 

p COS (p 

p sin cp 

0 

sin 2 ?) 

0 

U 

0 

cos cp/p 


0 

0 

U 

sin cp/p 


0 

o 

pa 6 cos cp 

pa 2 sin cp 

U 

(A15) 


The eigenvalues of this matrix are 


= sin I - r 2 sin ? 7 (Y 6*/X x yj + r| sin 2 ?) e 


(A16) 


where e represents one of the values U, U + a, or U - a. The value U is doubly 
degenerate. Equations (A8), (All), (A15), and (A16) can be combined to yield the fol- 
lowing expression for the eigenvalues of the matrix G: 

g = 1 - 4 k|( 1 - cos |) 2 + (1 - cos ?)) 2 ] - j A 2 e 2 - i 1(2 + cos \ + cos ?)) Ae (A17) 
The magnitudes of the eigenvalues g are given by the equation 
| g | 2 = (l . | A 2 e 2 ) 2 + iV 2 + cos 1 + cos 7) 2 A 2 e 2 

+ 4k|( 1 - cos i) 2 + (1 - cos i))^J ^4fc[(l - cos i) 2 + (1 - cos ?))^] - 2 + e 2 A 2 | (A18) 


Equation (A18) can be used to determine the minimum upper bound for At which 
satisfies the von Neumann condition for stability. In order that the damping terms (the 
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terms which are proportional to the damping coefficient k) will not increase the magni- 
tude of g, it is necessary that k not be negative and that the following inequality holds: 

4kJ(1 - cos l) 2 + (1 - cosrj)^j - 2 + e^A 2 < 0 

When this inequality is combined with the Courant-Friedrichs Lewy stability condition 
( e 2A 2 S l), an upper bound of l/32 is obtained for k. If the damping coefficient k 
satisfies the inequalities 

0 = k = 1/32 

equation (A 17) can be replaced by the inequality 

|g| 2 S 1 + (cos 1 + cos rj) 2 + (cos | + cos rj) - 3jA 2 e 2 

+ i A^e^ + 4 k j(l - cos l) 2 + (1 - cos T))^J(A^e 2 - l) (A19) 

The inequalities 

(cos | + cos r)) 2 = 2 (cos^f + cos 
cos I + cos rj = 2 
sin 2 | + sin 2 ?) § 2 


can be used to establish an upper bound for the right-hand side of inequality (A19). The 
resulting inequality can be written as 


1*1 


2 S 1 + i ^A^e 2 + 16 k J( 1 - cos l) 2 + (1 - cos j^e 2 - i(sin 2 f + sin 2 ?)) (A20) 


Since the quantity (- 2 e 2 + 16 k |(1 - cos 4)* + (1 - cos 7 ))^] J is positive definite, the 
magnitude of g is less than one only when the following inequality is satisfied: 


- O 

COS 4) + (1 - COS 




A 2 e 2 S A (sin 2 4 + sin 2 ?)) 


(A21) 
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From equation (A 16) it follows that 


e2s 


where 


'max 


Tjlsin §| + 


= u| + a 


r 2 |sinTj|Y|6'| 


+ r 2 |sin i)' 2 


'max 


(A22) 


The quantity cos p is defined so that 


cos p = 


sin | 


^sin 2 | + sin^ 


(A23) 


The inequalities (A21) and (A22) and equation (A23) can be combined to yield the 
inequality 


2a2<?2 5 2A 2 


sin 2 £ + sin 2 ?) 


r„ sin p y|5 t , „ „ 

rj cos p + — — J + r| sin*p 


e 2 

max 


i + 


( Y 6 ’Ax ) 2 + ( Y 6 ’Ax ) 2 ] 2 - 4r ? 


r 2 >e 2 £1 

r 2 r e max _ 


(A24) 


With this inequality and the definitions of A, r v r 2 , and U given in equations (A6) 
and (A14), it can be shown that At must satisfy the following condition if the 
von Neumann stability criterion is to be observed: 


At 5 


\j_ (* x AX) 2 + (6 AY) 2 / L J 


1 + 


(Y5* AX) 2 
AX^ 2 + (6 AY) 2 + ‘ 


l , (Y6’ AX) 2 

2 

2(\ x AX) (6 AY)" 

^A x AXj 2 + (6 AY) 2 


(\ x AX^ 2 + (6 AY) 2 


(A25) 


The smallest value of the right-hand side of this inequality which occurs anywhere in the 
flow field at a given time step constitutes the least upper bound for At for that partic- 
ular time step. It can be seen that the magnitude of the right-hand side of inequality (A25) 
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diminishes when the shock slope 6' becomes large (in other words, when the 
X,Y coordinate system becomes skewed). However, there is a compensating effect 
since the scale factor \ x and the shock-layer thickness 6 are generally large in 
regions where the coordinate system is skewed. 

For the calculations which are presented in this paper, it has been found that the 
least upper bound is located generally in regions where the scale factor A x or the 
shock -layer thickness 6 are small rather than in regions where the shock slope 5' 
is large. Since only convex bodies are treated in this paper, the smallest value of \ x 
occurs on the body surface at the segment with the smallest radius of curvature. In the 
case of bodies with sharp corners, the smallest pertinent value of \ x occurs on the 
line Y = AY where it subtends the corner. In general, the shock-layer thickness 5 
either attains or approaches its smallest value at the axis of symmetry. Therefore, the 
mesh spacings X x AX and 5 AY either approach or attain their smallest values on or 
near the stagnation streamline. 

The right-hand side of inequality (A25) can be simplified considerably when it is 
applied on the stagnation streamline since 6' = 0 on the axis of symmetry and Y = 0 
on the body surface. The expression for the upper bound of At which can be used on 
the stagnation streamline is 


At S 


min(A x AX, 6 AY) 
\J2 \/u 2 + (v - Y<5) 2 + 


(A 2 6) 


Since the shock speed is generally much smaller than the magnitude of the velocity, the 
term in inequality (A26) which is proportional to 6 can be neglected for practical pur- 
poses. Thus, the least upper bound for At is determined from the smallest value of 
the mesh spacings A x AX and 6 AY and the largest value of the sum V + a, where 
V is the total magnitude of the velocity. By using Bernoulli's equation, this sum can 
be written as 


V + a = a + J2 H - 


y - 1; 


(A 2 7) 


It can be shown that the maximum value of V + a for a given value of the total 
enthalpy H is 


(V + a} max = \Ar+DH = ^V raax 


(A28) 
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APPENDIX A - Concluded 


Since the total enthalpy does not change much during the course of a time-dependent 
calculation, the value of V max in the undisturbed free stream can be used safely in 
equation (A28). From the inequality (A26) and equation (A28), the following expression 
is obtained for the least upper bound of At: 


At S 


\]y + 1 V 


min(A x AX, 5 AY) 


max,°° 


(A29) 


The effect of the damping terms in equation (A17) is to prevent the finite -difference 
equations from becoming neutrally stable when one or more of the eigenvalues e vanish. 
Neutral stability occurs when the eigenvalue g has a magnitude of 1 ( j g | = l) . If the 
finite -difference equations are neutrally stable in terms of a linear analysis, it is quite 
possible that the nonlinear effects which have been neglected in the linear analysis will 
be sufficient to render the calculation unstable. It is seen from equations (A 14) and 
(A 16) that an eigenvalue e can vanish whenever one of the components u or (v - Y<5) 
vanishes or becomes sonic. 

It can be seen from equation (A17) that maximum damping occurs for L x = 2 AX 
and Ly = 2 AY although solutions with wavelengths of two mesh spacings are not 
observable. It is also seen that the damping decreases monotonically as the wavelengths 
of the disturbances increase. Therefore, the effect of the damping terms is to decrease 
the amplitude of short waves and to have very little influence on long waves. This is the 
desired effect since the physical solution being sought is generally characterized by 
wavelengths much longer than the mesh spacings being used, and solutions with wave- 
lengths of only a few mesh spacings are generally unwanted disturbances. 
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APPENDIX B 


PRANDT L -MEYER SOLUTION FOR THE TRANSIENT FLOW 
AT A SHARP CORNER 


In this appendix, the solution for compressible, time -dependent flow at a sharp 
corner is presented. It is shown that this solution is multiple valued and satisfies the 
Prandtl-Meyer equations. It should be emphasized that the solution is applicable only 
at the corner and is not valid in general in any neighborhood of the corner. 

When a perfect gas is being treated, equations (4) can be written as 


9p + iL + v + £- iiH + 
Bt A x Bx By \ x Bx 



+ jPu 


cos 8 


„ , ^r . sin 8 
+ pV ‘Aj + J -Y 


= 0 




9H + iL9H + v^ + A.^ + ^uv = 0 

at a x ax By pa x ax a x 


By [ u 9v 9v 

at Ax 9x 9 y 


+ 19E-^1u 2 = 0 

P By A x 


9E + JL®E + v 5E-2^f5fi + -!i-^ + v^ = 0 

at A x ax ay P yat a x ax Byj 


> (Bl) 


It should be noted that the curvature of that part of the reference line which subtends the 
corner is 



Let x c max and x c m ^ n be the maximum and minimum values, respectively, of the 
coordinate x associated with the corner. If the angle 8 for x £ m ^ n is denoted 
by 0 C , then the function 9(x) can be written as 


8 = 


9 c + 


x c,min 

yb 


for x c m i n = x = x c max - Let a new independent variable R be defined such that 


R = y + y b 
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This coordinate is measured from the corner along lines perpendicular to the reference 
surface. It should be noted that X x and r are written in terms of R and 0 as 


X 


x 


R_ 


r = r^ + R sin 9 


where r b is the perpendicular distance from the corner to the axis of symmetry. 
Equations (Bl) are written in terms of R, 9, and t as 


9P + y_9P + v^ + £.9H + p 9v +j p u cosj 

at R d9 9R R 99 3R r b + R sin 9 


pvfl + j *gL» ' \ = 0 

\R J r b + R sin 9 J 


9H + y.9H + v^. + -L9E + HZ =0 

at R a e 9R pR ae r 


2l + H.9Z + v ^ + I^ P-^ =0 


at R 99 


3R p3R R 


> (B2) 


+ u ^2 + v ®E - + H.M + v^ = 0 

at R 30 3R P \at R 99 9R) 


J 


Let the flow properties p, p, u 5 and v in the vicinity of the corner be represented by 
functions of the form 


F(t,R, 0) = F c (t,0) + R a(t ’ 0) F R (t,R,0) 


(B3) 


where F R (t,R,0) is bounded for small values of R and the exponent is greater than or 
equal to zero. If these functions are substituted into equations (B2), and if these equa- 
tions are then multiplied by R and the limit is taken as R approaches zero, the terms 
involving a and F R are eliminated and the following equations are obtained: 



(B4a) 


(B4b) 
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u l^- 

c c\88 


3p c 

p c'iF‘ yp 


u c) = 0 


c 9 e 


(B4c) 

(B4d) 


These equations are the Prandtl-Meyer equations where, by assumption, p c , P c , u c , 
and v c are functions of the time t as well as of the angle 9. 

The solution for the flow properties p, p, u, and v at the corner are of the 

form 


F(t,x,y) = F c (t, 9) 


where x ^ S x S x air and y = -y b . If x* £ x c min is the smallest value of x 
associated with the corner for which u = a, if 9* is the angle 8 for x = x*, and if 
the flow properties and the angle 9 for x = x c min are denoted by p c , p c , u c , 
and 9 C , respectively, then the solution at the corner which satisfies equations (B4) can 
be written as 


P = P C 


P = ■ P, 


u = u^, cos 


( s - «c) 


v = U£ sin(0 - Oq} 


(B5) 


for x . = x = x* and 
c,min 


P = P C 


(COS 


y - l 

y + 1 


(8 - 9*) 


- ft ttK - l ) si " - e *> 


P = P c < cos 


r 

l 


y - l 

y + 1 


(9 - 9*) 




V-ll(9 - 9 *) 
y + 1 


_2y 
iY- 1 


_ 2 _ 

\y-i 


(B6a) 


(B6b) 
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APPENDIX B - Concluded 



(B6c) 

(B6d) 


for x„ ^ x £ x*. The quantities a n and M,-, in equations (B6) are the speed of 

C , IilcLX v/ U 

sound and the Mach number, respectively, at the corner for x = x c m ^ n . 
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Figure 5«- Surface pressure and velocity distributions for flat-face cylinder -with r c/ r t> : 










0 .2 .4 .6 .8 1.0 0 .2 .4 .6 .8 1.0 

Y Y 


(a) Pressure. (b) Density. 

Figure 6.- Present results for distributions of flow properties across shock layer 
for flat -face cylinder with r^r^ = 0.4. y = 1.4; ^ = 2.49- The coordinate Y 

has values of 1 and 0 at the shock wave and body surface, respectively. 
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(c) Tangential velocity component. (d) Normal velocity component. 

Figure 6.- Concluded. 
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(a) Pressure. (t>) Density. 

Figure 9»- Present results for distributions of flow properties across shock layer 
for flat-face cylinder with sharp comer, 7 = 1. b; ^ = 2.81. The coordinate Y 
has values of 1 and 0 at the shock wave and "body surface, respectively. 
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Effect of parabolic free-stream velocity distribution on shock- 
sonic-line locations for spherically blunted 60 ^ cone with 
>.25. y = 1.4; 14 ^ = 10 . 
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